Fort Bend County District Boundaries

All shapefiles were projected into the Texas mapping projection for the purposes of creating maps and calculating the physical properties of the districts’ shape. The offices we are focusing on are county commissioner, justice of the peace, and constable. The geographic boundaries for these offices are referred to as districts or precincts. Some counties will refer to them as commissioner precinct 1, 2, etc. as opposed to voting precincts. Throughout the report, district and precinct may be used interchangeably. Change this to be consistent with how the county refers to it.

######### GEOPROCESSING

#### PRECINCTS

## EXTRACT PRECINCT & DISTRICT IDs FROM COUNTY
county_precs_dists <- precinct_district_cw %>%
  filter(county_fips == cnty_fips) %>%
  select(van_precinct_name,
         dnc_precinct_id,
         "Commish" = `County Commissioner`,
         JP)

## CLEAN GEOMETRY & JOIN TO PRECINCTS FILE
county_precincts_prj <- tx_precincts %>%
  filter(CNTY == cnty_fips) %>%
  select(PREC, dnc_id, c_name) %>%  
  st_make_valid(.) %>%
  st_buffer(., dist = 0) %>%
  left_join(county_precs_dists, 
            by = c("dnc_id" = "dnc_precinct_id"))


#### CREATE COMMISSIONERS DISTRICTS
cnty_commish_dists_prj <- county_precincts_prj %>%
  select(district_type = "Commish") %>% 
  group_by(district_type) %>%
  summarize(n_precincts = n()) %>%
  mutate(area = st_area(.),
         perimeter = st_perimeter(.),
         district_type = paste0("CC-", district_type)) %>%
  st_buffer(., 0.5) %>%
  st_difference(., .) %>%
  group_by(district_type, n_precincts, area, perimeter) %>%
  summarize(geometry = st_union(geometry))

  # check for additional intersections & test diff buffer widths for cleaning
#ggplot() + geom_sf(data=cnty_commish_dists_prj, aes(fill = district_type))

#### CREATE JP DISTRICTS
cnty_JP_dists_prj <- county_precincts_prj %>%
  select(district_type = "JP") %>% 
  group_by(district_type) %>% 
  summarize(n_precincts = n()) %>%
  mutate(area = st_area(.),
         perimeter = st_perimeter(.),
         district_type = paste0("JP-", district_type)) %>%
  st_buffer(., 0.5) %>%
  st_difference(., .) %>%
  group_by(district_type, n_precincts, area, perimeter) %>%
  summarize(geometry = st_union(geometry))
  
  # check for additional intersections & test diff buffer widths for cleaning
#ggplot() + geom_sf(data=cnty_JP_dists_prj, aes(fill = district_type))

Some text about Map 1 about here.

Map 1: Reference Map of Commissioners Districts

Some text about Map 2 about here.

Map 2: Reference Map of JP & Constable Districts

## GEOPROCESSING OF OTHER SPATIAL FEATURES FOR MAPS 
    #If you need other spatial features for the maps geoprocess in this chunk.

Map 3 is an interactive map of the county commissioner districts/precincts in Fort Bend County. Click on the zoom in button to see what physical features and landmarks are within each district. One voting precinct is missing a county commissioner assignment in the TDP all-geography-overlap.csv. This voting precinct is labeled as NA in the map.

## COMMISH INTERACTIVE MAPPING
commish_pal <- colorFactor(
  palette = "Spectral",
  domain = cnty_commish_dists_prj$district_type
)

leaflet() %>%
  addTiles() %>% 
  addPolygons(
    data = cnty_commish_dists_prj %>% st_transform(prj_webmap),
    fillColor = ~commish_pal(district_type),
    color = "white",
    opacity = 0.5,
    fillOpacity = 0.5,
    weight = 1) %>% 
  addLegend(
    data = cnty_commish_dists_prj,
    pal = commish_pal,
    values = ~district_type,
    title = "Map 3: County Commissioner District<br>Interactive Map"
  )

Map 4 is an interactive map of the JP and constables districts/precincts in Fort Bend County. Click on the zoom in button to see what physical features and landmarks are within each district. You’ll notice that the boundaries of county commissioner districts are identical to the JP/Constable precincts. There are four positions for all three types of offices in the county. As in Map 3 with the commissioner’s district, one voting precinct is missing a JP assignment in the TDP all-geography-overlap.csv. This precinct is labeled as NA in the map.

## JP/CONSTABLE INTERACTIVE MAPPING
JP_pal <- colorFactor(
  palette = "Accent",
  domain = cnty_JP_dists_prj$district_type
)

leaflet() %>%
  addTiles() %>% 
  addPolygons(
    data = cnty_JP_dists_prj %>% st_transform(prj_webmap),
    fillColor = ~JP_pal(district_type),
    color = "white",
    opacity = 0.5,
    fillOpacity = 0.5,
    weight = 1) %>% 
  addLegend(
    data = cnty_JP_dists_prj,
    pal = JP_pal,
    values = ~district_type,
    title = "Map 4: JP/Constable District<br>Interactive Map"
  )

Properties of the District’s Boundaries

The area and perimeter of each district/precinct was calculated. Then, several compactness measures were computed, which are displayed below in Table 1. For all four measures, a score closer to 1 indicates higher compactness.

Note: County commissioners and JP/Constable districts/precinct have the same boundaries. Also, be mindful of the NA district/precinct at this time.

Among the county commissioners’ districts, district 1 stands out as the least compact (i.e., most irregular). District 1 scores higher on convexity than all of the other measures, which indicates that the district has a convex shape. However, the elongated maximum width of the district and its area set it apart as particularly irregular.

Similarly, the JP districts tend to be more convex in shape. However, a handful score high on elongation, like JP districts 2, 5, and 8.

##### CALCULATE COMPACTNESS MEASURES
polsby_popper = function(poly1) {
  return(drop_units(
    4 * pi * st_area(poly1) / st_length(st_boundary(poly1))^2))
}

schwartzberg = function(poly1) {
  return(polsby_popper(poly1)^-0.5)
}

inv_schwartzberg = function(poly1) {
  return(polsby_popper(poly1)^0.5)
}

area_convex_hull = function(poly1, ch = NULL) {
  if (is.null(ch)) {
    ch = st_convex_hull(st_geometry(poly1))
  }
  return(drop_units(st_area(poly1) / st_area(ch)))
}

reock = function(poly1, mbc = NULL) {
  if (is.null(mbc)) {
    mbc = st_minimum_bounding_circle(st_convex_hull(st_geometry(poly1)))
  }
  return(drop_units(st_area(poly1) / st_area(mbc)))
}

cnty_districts_land_table <- rbind(cnty_commish_dists_prj, cnty_JP_dists_prj) %>%
  mutate(polsby_popper = polsby_popper(geometry),
         schwartzberg = schwartzberg(geometry),
         inv_schwartzberg = inv_schwartzberg(geometry),
         reock = reock(geometry),
         area_convex_hull = area_convex_hull(geometry))

knitr::kable(cnty_districts_land_table %>%
               st_set_geometry(NULL) %>%
               select(-schwartzberg) %>%
               rename(`District Type` = "district_type",
                       Precincts = n_precincts,
                      `polsby-popper` = "polsby_popper",
                      `area convex hull` = "area_convex_hull",
                      `inverse schwartzberg` = "inv_schwartzberg"
             ), caption = "Table 1: Compactness Scores by District")
Table 1: Compactness Scores by District
District Type Precincts area perimeter polsby-popper inverse schwartzberg reock area convex hull
CC-1 39 1540260962.7 [m^2] 280271.357 [m] 0.2464521 0.4964394 0.4605032 0.7737425
CC-2 45 131088346.8 [m^2] 108477.520 [m] 0.1401046 0.3743055 0.1787794 0.4850524
CC-3 33 434979532.1 [m^2] 196115.840 [m] 0.1421729 0.3770582 0.2524688 0.6589315
CC-4 42 178549550.1 [m^2] 104981.224 [m] 0.2036752 0.4513039 0.2820256 0.7013570
CC-NA 1 298291.6 [m^2] 2467.829 [m] 0.6164694 0.7851557 0.4219893 0.8633141
JP-1 39 1540260962.7 [m^2] 280271.357 [m] 0.2464521 0.4964394 0.4605032 0.7737425
JP-2 45 131088346.8 [m^2] 108477.520 [m] 0.1401046 0.3743055 0.1787794 0.4850524
JP-3 33 434979532.1 [m^2] 196115.840 [m] 0.1421729 0.3770582 0.2524688 0.6589315
JP-4 42 178549550.1 [m^2] 104981.224 [m] 0.2036752 0.4513039 0.2820256 0.7013570
JP-NA 1 298291.6 [m^2] 2467.829 [m] 0.6164694 0.7851557 0.4219893 0.8633141

Population Characteristics

Table 2, shown below, summarizes the growth in Fort Bend County from the last decennial Census in 2010 through 2019. Fort Bend has grown substantially over the decade. The county has added just over 180,000 people, which is a 30% change from 2010. The citizen-voting age population has grown by 13% over this period. Be mindful that the 2019 data are estimates, so the CVAP population counts may be an underestimate. This value should be validated using other datasets.

## 2010
census_2010 <- readRDS("census_2010.rds")
census_2010_county <- census_2010 %>%
  mutate(county_fips = as.numeric(substr(cbg, 3, 5))) %>%
  filter(county_fips == cnty_fips) %>%
  summarize(total_population_2010 = sum(total_population, na.rm =T),
            cvap_2010 = sum(cvap, na.rm =T)) %>%
  pivot_longer(ends_with("_2010"),
    names_to = "variable",
    values_to = "Total, 2010") %>% 
  mutate(variable = str_remove(variable, "_2010"),
         variable = str_replace_all(variable, "_", " "))

## 2019
acs19_full <- readRDS("acs19_full.rds")
acs19_full_county <- acs19_full %>%
  mutate(county_fips = as.numeric(substr(cbg, 3, 5))) %>%
  filter(county_fips == cnty_fips) 
  

## Time Table
knitr::kable(acs19_full_county %>%
  summarize(total_population_2019 = sum(pop, na.rm =T),
            cvap_2019 = sum(CVAP, na.rm =T)) %>%
  pivot_longer(ends_with("_2019"),
    names_to = "variable",
    values_to = "Total, 2019") %>% 
  mutate(variable = str_remove(variable, "_2019"),
         variable = str_replace_all(variable, "_", " ")) %>%
  left_join(census_2010_county, by = c("variable")) %>%
  relocate(variable, `Total, 2010`, `Total, 2019`) %>%
  mutate(`Total Change` = `Total, 2019` - `Total, 2010`,
         `% Change` = `Total Change` / `Total, 2010`,
         `% Change` =  percent(`% Change`,
                               accuracy = 0.01,
                               scale = 100,
                               label = "%"),
         `Total Change` = format(`Total Change`, big.mark   = ","),
         `Total, 2010` = format(`Total, 2010`, big.mark   = ","),
         `Total, 2019` = format(`Total, 2019`, big.mark = ",")), 
  caption = "Table 2: Fort Bend County Population 2010 & 2019")
Table 2: Fort Bend County Population 2010 & 2019
variable Total, 2010 Total, 2019 Total Change % Change
total population 585,375 765,394 180,019 30.75%
cvap 411,540 466,249 54,709 13.29%
## CVAP 2019 Summary
CVAP19_summary <- acs19_full_county %>%
  select(c(3:19)) %>%
  select(c(-6:-7, -9:-12)) %>%
  summarise_if(is.numeric, sum, na.rm = T) %>%
  relocate(`High school graduate (includes equivalency)`,
           .after = `No high school diploma`) %>%
  pivot_longer(cols = c(1:11),
               names_to = "Group", 
               values_to = "Total") %>%
  mutate(`% of CVAP` = percent(Total / first(Total),
                           accuracy = 0.01,
                           scale = 100,
                           label = "%" ),
         Total = format(Total, big.mark = ","))

knitr::kable(CVAP19_summary,
            caption = "Table 2.1: Summary of Citizen Voting-Age Demographics")
Table 2.1: Summary of Citizen Voting-Age Demographics
Group Total % of CVAP
CVAP 466,249 100.00%
18 to 29 years 91,881 19.71%
30 to 44 years 130,043 27.89%
45 to 64 years 171,000 36.68%
65 years and over 73,325 15.73%
Income in the past 12 months below poverty level 25,960 5.57%
Income in the past 12 months at or above the poverty level 435,787 93.47%
No high school diploma 36,221 7.77%
High school graduate (includes equivalency) 89,992 19.30%
Some college or Associates 143,136 30.70%
Bachelor’s or Graduate degree 196,900 42.23%
## Race 2019 Summary
Race19_summary <- acs19_full_county %>%
  select(c(2, 20:27)) %>%
  rename(`Total population` = "pop") %>%
  summarise_if(is.numeric, sum, na.rm = T) %>%
  pivot_longer(cols = c(1:9),
               names_to = "Group", 
               values_to = "Total") %>%
  mutate(`% of Population` = percent(Total / first(Total),
                           accuracy = 0.01,
                           scale = 100,
                           label = "%" ),
         Total = format(Total, big.mark = ","),
         Group = str_replace(Group, "_", " "))

knitr::kable(Race19_summary,
            caption = "Table 2.3: Summary of Racial Demographics")
Table 2.3: Summary of Racial Demographics
Group Total % of Population
Total population 765,394 100.00%
White Not Latinx 253,263 33.09%
Black Not Latinx 153,972 20.12%
Two or More Not Latinx 1,474 0.19%
Other Not Latinx 1,559 0.20%
NHPI Not Latinx 396 0.05%
Asian Not Latinx 153,245 20.02%
American Indian Not Latinx 1,713 0.22%
Any race Latinx 187,500 24.50%
## Household 2019 Summary
HH19_summary <- acs19_full_county %>%
  select(c(28:45)) %>%
  relocate(limited_english, .after = Households) %>%
  rename(`Limited English Proficiency` = "limited_english") %>%
  summarise_if(is.numeric, sum, na.rm = T) %>%
  pivot_longer(cols = c(1:18),
               names_to = "Group", 
               values_to = "Total") %>%
  mutate(`% of Population` = percent(Total / first(Total),
                           accuracy = 0.01,
                           scale = 100,
                           label = "%" ),
         Total = format(Total, big.mark = ","),
         Group = str_replace(Group, "_", " "))

knitr::kable(HH19_summary,
            caption = "Table 2.4: Summary of Household Demographics")
Table 2.4: Summary of Household Demographics
Group Total % of Population
Households 237,883 100.00%
Limited English Proficiency 15,073 6.34%
Less than $10,000 7,437 3.13%
$10,000 to $14,999 4,244 1.78%
$15,000 to $19,999 4,766 2.00%
$20,000 to $24,999 6,146 2.58%
$25,000 to $29,999 6,185 2.60%
$30,000 to $34,999 6,547 2.75%
$35,000 to $39,999 6,890 2.90%
$40,000 to $44,999 7,148 3.00%
$45,000 to $49,999 6,246 2.63%
$50,000 to $59,999 15,132 6.36%
$60,000 to $74,999 20,550 8.64%
$75,000 to $99,999 29,748 12.51%
$100,000 to $124,999 27,914 11.73%
$125,000 to $149,999 20,743 8.72%
$150,000 to $199,999 29,417 12.37%
$200,000 or more 38,770 16.30%
#### GEOPROCESSING OF BLOCK GROUPS
cnty_bgs_prj <- county_bgs %>%
  dplyr::select(GEOID) %>%
  st_transform(st_crs(cnty_commish_dists_prj)) %>%
  left_join(acs19_full_county, by = c("GEOID" = "cbg")) %>%
  dplyr::select(-GEOID, -county_fips)

  
#### APPORTION DEMOGRAPHIC POPULATIONS TO DISTRICTS
cnty_commish_apportion <- st_interpolate_aw(cnty_bgs_prj,
                                           cnty_commish_dists_prj,
                                           extensive = T)

cnty_JP_apportion <- st_interpolate_aw(cnty_bgs_prj,
                                       cnty_JP_dists_prj,
                                       extensive = T)

In Table 2.5, the NA district/precinct shows up as #5.

Table 2.5: County Commissioner’s CVAP Apportionment

  # CVAP apportionment - CC
cnty_commish_apportion_cvap <- cnty_commish_apportion %>%
  dplyr::select(Group.1, pop:X65.years.and.over,
                High.school.graduate..includes.equivalency.,
                Income.in.the.past.12.months.below.poverty.level:Bachelor.s.or.Graduate.degree) %>%
  relocate(No.high.school.diploma, .after = X65.years.and.over) %>%
  relocate(Some.college.or.Associates, .after = High.school.graduate..includes.equivalency.) %>%
  relocate(Bachelor.s.or.Graduate.degree, .after=Some.college.or.Associates) %>%
  rename(District = "Group.1",
         total_population = "pop") %>%
  st_set_geometry(NULL) %>%
  pivot_longer(cols = c(2:13),
               names_to = "Group",
               values_to = "estimated_count") %>%
  mutate(Group = str_remove(Group, "X"),
         Group = str_replace_all(Group, "/.", " "),
         District = paste0("CC-0", District)) %>%
  pivot_wider(id_cols = District,
              names_from = "Group",
              values_from = "estimated_count") %>%
  mutate_each(funs(prettyNum(.,  big.mark=",")))

names(cnty_commish_apportion_cvap) <-  gsub("[[:punct:]]",
                                            " ",
                                            names(cnty_commish_apportion_cvap))

rmarkdown::paged_table(cnty_commish_apportion_cvap)

In Table 2.5, the NA district/precinct shows up as #5.

Table 2.6: Apportionment of Race by Commissioner’s District

 # race apportionment - CC
cnty_commish_apportion_race <- cnty_commish_apportion %>%
  dplyr::select(Group.1, pop, White.Not_Latinx:Any.race.Latinx) %>%
  rename(District = "Group.1",
         total_population = "pop") %>%
  st_set_geometry(NULL) %>%
  pivot_longer(cols = c(total_population:Any.race.Latinx),
               names_to = "Group",
               values_to = "estimated_count") %>%
  mutate(Group = str_remove(Group, "X"),
         Group = str_replace_all(Group, "/.", " "),
         District = paste0("CC-0", District)) %>%
  pivot_wider(id_cols = District,
              names_from = "Group",
              values_from = "estimated_count") %>%
  mutate_each(funs(prettyNum(.,  big.mark=",")))  

names(cnty_commish_apportion_race) <-  gsub("[[:punct:]]",
                                            " ",
                                            names(cnty_commish_apportion_race))


rmarkdown::paged_table(cnty_commish_apportion_race,list(rows.print = 4))

Table 2.7: JP/Constable CVAP Apportionment

  # CVAP apportionment - JP
cnty_JP_apportion_cvap <- cnty_JP_apportion %>%
  dplyr::select(Group.1, pop:X65.years.and.over,
                High.school.graduate..includes.equivalency.,
                Income.in.the.past.12.months.below.poverty.level:Bachelor.s.or.Graduate.degree) %>%
  relocate(No.high.school.diploma, .after = X65.years.and.over) %>%
  relocate(Some.college.or.Associates, .after = High.school.graduate..includes.equivalency.) %>%
  relocate(Bachelor.s.or.Graduate.degree, .after=Some.college.or.Associates) %>%
  rename(District = "Group.1",
         total_population = "pop") %>%
  st_set_geometry(NULL) %>%
  pivot_longer(cols = c(2:13),
               names_to = "Group",
               values_to = "estimated_count") %>%
  mutate(Group = str_remove(Group, "X"),
         Group = str_replace_all(Group, "/.", " "),
         District = paste0("JP-0", District)) %>%
  pivot_wider(id_cols = District,
              names_from = "Group",
              values_from = "estimated_count") %>%
  mutate_each(funs(prettyNum(.,  big.mark=",")))  

names(cnty_JP_apportion_cvap) <-  gsub("[[:punct:]]",
                                       " ",
                                       names(cnty_JP_apportion_cvap))

rmarkdown::paged_table(cnty_JP_apportion_cvap, list(rows.print = 8))

Table 2.8: Apportionment of Race by JP/Constable District

 # race apportionment - JP
cnty_JP_apportion_race <- cnty_JP_apportion %>%
  dplyr::select(Group.1, pop, White.Not_Latinx:Any.race.Latinx) %>%
  rename(District = "Group.1",
         total_population = "pop") %>%
  st_set_geometry(NULL) %>%
  pivot_longer(cols = c(total_population:Any.race.Latinx),
               names_to = "Group",
               values_to = "estimated_count") %>%
  mutate(Group = str_remove(Group, "X"),
         Group = str_replace_all(Group, "/.", " "),
         District = paste0("JP-0", District)) %>%
  pivot_wider(id_cols = District,
              names_from = "Group",
              values_from = "estimated_count") %>%
  mutate_each(funs(prettyNum(.,  big.mark=",")))  

names(cnty_JP_apportion_race) <-  gsub("[[:punct:]]",
                                        " ",
                                        names(cnty_JP_apportion_race))

rmarkdown::paged_table(cnty_JP_apportion_race, list(rows.print = 8, cols.print = 9))

2020 Election Results & Political Partisanship

Election Results 2020 Need to obtain election results from county website.

Partisanship Summary Pull data from partisanship model.

Redistricting - Demo of Algorithm for Creating New Districts

# CREATE ROOK NEIGHBORS
precinct_adjacency <- st_rook(county_precincts_prj) # this IDs adjacent precincts
  # use rooks, not queen neighbors for contiguity

# iNTERPOLATE BLOCK GROUPS TO PRECINCTS
precincts_demography <- cnty_bgs_prj %>%
  select(pop, CVAP, `White Not_Latinx`:`Any race Latinx`) %>%
  st_interpolate_aw(., 
                    county_precincts_prj,
                    extensive = T) %>%
  st_set_geometry(NULL)

# DATA FRAME FOR REDISTRICTING JP/CONSTABLE
precincts_JPConst_data <- county_precincts_prj %>%
  select(-Commish) %>%
  st_set_geometry(NULL) %>%
  cbind(., precincts_demography) %>%
  mutate(PREC = as.numeric(PREC),
         JP = as.numeric(JP),
         JP = replace_na(JP, 3))
# JOIN IN THE POLITICAL OUTCOMES DATA HERE
# %>%
#   left_join(., constable_2020_all, by = c("PREC" = "precinct")) %>%
#   left_join(., partisanship, by = c("dnc_precinct_id")) %>%
#   select(-PREC, -Group.1, -democrat:-republican) %>%
#   mutate(dems = highd + lowd,
#          repubs = total_reg_voters - dems)
# Redistricting with tpop r1
alg_mcmc_JPConst_initcd <- redist.mcmc(adjobj = precinct_adjacency,
                                       popvec = precincts_JPConst_data$pop, 
                                       initcds = precincts_JPConst_data$JP, 
                                       nsims = 10000,
                                       popcons = 0.10,
                                       maxiterrsg = 45000, 
                                       savename = paste0(county, 
                                                         "_alg_mcmc_JPConst_redist")
                            )
## 
## ==================== 
## redist.mcmc(): Automated Redistricting Simulation Using
##          Markov Chain Monte Carlo
## 
## Preprocessing data.
## 
## 10 percent done.
## Metropolis acceptance ratio: 0.968969
## 
## 20 percent done.
## Metropolis acceptance ratio: 0.967484
## 
## 30 percent done.
## Metropolis acceptance ratio: 0.968656
## 
## 40 percent done.
## Metropolis acceptance ratio: 0.964991
## 
## 50 percent done.
## Metropolis acceptance ratio: 0.963793
## 
## 60 percent done.
## Metropolis acceptance ratio: 0.961994
## 
## 70 percent done.
## Metropolis acceptance ratio: 0.962566
## 
## 80 percent done.
## Metropolis acceptance ratio: 0.962245
## 
## 90 percent done.
## Metropolis acceptance ratio: 0.960885
## 
## 100 percent done.
## Metropolis acceptance ratio: 0.960396

It looks like most of the plans developed by the algorithm are within the 10% population deviation requirement.

summary(alg_mcmc_JPConst_initcd$distance_parity)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.004436 0.058420 0.074422 0.071836 0.087307 0.113819
  # At this point you'd want to extract all those under 0.1 as potential plans.

If we compare the plan that does the best job at equalizing population across the four districts to the current plan, you’ll see that only a few precincts were moved around from district 2 to 3 in particular. n.b. 0 is technically district 1, 1 is 2, etc. The redistricting algorithm prefers to 0-index each unit.

# PRECINCTS BY NEW DISTRICTS
table(alg_mcmc_JPConst_initcd$partitions[, 1])
## 
##  0  1  2  3 
## 39 48 31 42
# PRECINCTS IN CURRENT DISTRICTS
table(precincts_JPConst_data$JP)
## 
##  1  2  3  4 
## 39 45 34 42
# CONVERT TO DATAFRAME AND MAP 
JPConst_mcmc_init_partitions <- as.data.frame(alg_mcmc_JPConst_initcd$partitions[,1])
JPConst_mcmc_partitions_init_comparison <- JPConst_mcmc_init_partitions %>%
  rename(partitions = "alg_mcmc_JPConst_initcd$partitions[, 1]") %>%
  mutate(partitions = as.factor(partitions + 1)) %>% 
  cbind(county_precincts_prj) %>%
  st_as_sf(.)

current_init_JPConst <- ggplot() + 
  geom_sf(data=JPConst_mcmc_partitions_init_comparison, aes(fill = as.factor(JP))) +
  theme(axis.title.x=element_blank(),
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.title.y=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks.y=element_blank())

new_init_JPConst <- 
  ggplot() + geom_sf(data=JPConst_mcmc_partitions_init_comparison, aes(fill = partitions)) + 
  theme(axis.title.x=element_blank(),
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.title.y=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks.y=element_blank())

ggarrange(current_init_JPConst, new_init_JPConst, 
          labels = c("Current Boundaries", "New Boundaries"),
          ncol = 2, nrow = 1)

After testing out some different plans, you could also run some segregation analyses on subgroups to see how the different plans influence the distribution of subgroups, e.g., Black population, AAPI population, democratic voters, limited ELP, etc. This analysis uses the disimilarity index. These can indicate which precincts have more/less of one group than another.

In the segregation indices, the values indicate how much of the population in a precinct need to move to another neighborhood (precinct) in order to achieve integration in the precinct. For example, in Fort Bend, on the low end, 22% of Black residents would need to move to another precinct in order to achieve integration with non-Black residents. At the higher end, 26% of Black residents would need to move to another precinct to be distributed the same way as the non-Black population in the precinct. We could interpret the precincts with a 0.26 dissimilarity index as the precincts where Black residents are most separated from other racial groups.

# SEG ANALYSIS
seg_JPConst_initcd <- redist.segcalc(algout = alg_mcmc_JPConst_initcd, 
                      grouppop = precincts_demography$Black.Not_Latinx,
                      fullpop = precincts_demography$pop)

seg_JPConst_initcd_df <-as.data.frame(seg_JPConst_initcd[1:160])

JPConst_mcmc_partitions_init_comparison <- cbind(JPConst_mcmc_partitions_init_comparison,
                                                 seg_JPConst_initcd[1:160])
ggplot() + 
  geom_sf(data=JPConst_mcmc_partitions_init_comparison,
          aes(fill = seg_JPConst_initcd.1.160.))

You can also apply constraints to the algorithm. Some potential ideas for constraints include constraints on communities of interest with subpopulations like race groups, voting rights act criteria, compactness, and geographic features (e.g., municipal boundaries, rivers, highways).

## VRA & SIMILARITY ON BLACK POPULATION ACROSS DISTRICTS
alg_mcmc_JPConst_initcd_BlackVRA <- redist.mcmc(adjobj = precinct_adjacency,
                                       popvec = precincts_demography$pop,
                                       initcds = precincts_JPConst_data$JP,
                                       nsims = 10000,
                                       popcons = 0.10,
                                       grouppopvec = precincts_JPConst_data$Black.Not_Latinx,
                                       constraint = c("vra", "similarity"),
                                       constraintweights = c(5, 3), 
                                       maxiterrsg = 45000,
                                       savename = paste0(county, 
                                                         "_alg_mcmc_JPConst_redist_BlackVRA")
                                       )
## 
## ==================== 
## redist.mcmc(): Automated Redistricting Simulation Using
##          Markov Chain Monte Carlo
## 
## Preprocessing data.
## 
## 10 percent done.
## Metropolis acceptance ratio: 0.962963
## 
## 20 percent done.
## Metropolis acceptance ratio: 0.963482
## 
## 30 percent done.
## Metropolis acceptance ratio: 0.960987
## 
## 40 percent done.
## Metropolis acceptance ratio: 0.96049
## 
## 50 percent done.
## Metropolis acceptance ratio: 0.962192
## 
## 60 percent done.
## Metropolis acceptance ratio: 0.959493
## 
## 70 percent done.
## Metropolis acceptance ratio: 0.958994
## 
## 80 percent done.
## Metropolis acceptance ratio: 0.959495
## 
## 90 percent done.
## Metropolis acceptance ratio: 0.958329
## 
## 100 percent done.
## Metropolis acceptance ratio: 0.958196
# VIEW SUMMARY STATS OF NEW PLANS 
summary(alg_mcmc_JPConst_initcd_BlackVRA$distance_parity)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.001852 0.058064 0.074973 0.071817 0.087393 0.111624
table(alg_mcmc_JPConst_initcd_BlackVRA$partitions[, 3])
## 
##  0  1  2  3 
## 40 47 32 41
# CONVERT TO DATAFRAME AND MAP 
JPConst_mcmc_init_constraints_partitions <- as.data.frame(alg_mcmc_JPConst_initcd_BlackVRA$partitions[,2])
JPConst_mcmc_partitions_init_comparison <- JPConst_mcmc_init_constraints_partitions %>%
  rename(partitions = "alg_mcmc_JPConst_initcd_BlackVRA$partitions[, 2]") %>%
  mutate(partitions = as.factor(partitions + 1)) %>% 
  cbind(county_precincts_prj) %>%
  st_as_sf(.)

new_init_JPConst_constraints <- 
  ggplot() + geom_sf(data=JPConst_mcmc_partitions_init_comparison, aes(fill = partitions)) +   theme(axis.title.x=element_blank(),
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.title.y=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks.y=element_blank())


black_pop_choropleth <- cbind(county_precincts_prj, precincts_JPConst_data) %>%
  ggplot() + 
  geom_sf(aes(fill = Black.Not_Latinx))  +
  scale_fill_continuous(trans = 'reverse') + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

ggarrange(current_init_JPConst, new_init_JPConst_constraints, black_pop_choropleth,
          labels = c("Current Boundaries", "New Boundaries",
                     "Black Population Estimates"),
          ncol = 2, nrow = 2)